N90- 23027 

A Finite Element Based Method for 
Solution of Optimal Control Problems 

Robert R. Bless \ Dewey H. Hodges 1 2 , and Anthony J. Calise 3 

School of Aerospace Engineering 
Georgia Institute of Technology, Atlanta, GA 30332 

Abstract 

A temporal finite element based on a mixed form of the Hamiltonian weak principle is presented for optimal 
control problems. The mixed form of this principle contains both states and costates as primary variables that 
are expanded in terms of elemental values and simple shape functions. Unlike other variational approaches to 
optimal control problems, however, time derivatives of the states and costates do not appear in the governing 
variational equation. Instead, the only quantities whose time derivatives appear therein are virtual states and 
virtual costates. Also noteworthy among characteristics of the finite element formulation is the fact that in the 
algebraic equations which contain costates, they appear linearly. Thus, the remaining equations can be solved 
iteratively without initial guesses for the costates; this reduces the size of the problem by about a factor of two. 
Numerical results are presented herein for an elementary trajectory optimization problem which show very good 
agreement with the exact solution along with excellent computational efficiency and self-starting capability. The 
goal of this work is to evaluate the feasibility of this approach for real-time guidance applications. To this end, a 
simplified two-stage, four-state model for an advanced launch vehicle application is presented which is suitable 
for finite element solution. 


Introduction 

Future space transportation and deployment needs are critically dependent on the development of reliable 
and economical launch vehicles that will provide flexible, routine access to orbit. A particular requirement now 
receiving attention is that for an advanced technology heavy-lift vehicle. Future space transportation Systems 
will need to place large payloads - 100,000 to 150,000 pounds - into low Earth orbit at an order of magnitude 
lower cost per pound. Such systems will also require on-board algorithms that maximize system performance 
as measured by autonomy, mission flexibility, in-flight adaptability, reliability, accuracy and payload capability. 
They must be computationally efficient, robust, self-starting, and capable of functioning independently of ground 
control. Also, the algorithms must be designed with the anticipation that the launch vehicle will undergo 
evolutionary growth [1]. 

One approach to optimal guidance consists of repeatedly solving a two-point boundary- value problem that 
results from the traditional necessary conditions for optimality in an optimal control problem formulation. The 
vehicle state at discrete instants of time along a trajectory can be viewed as a new starting condition, and the 
remainder of the trajectory is reoptimized for that condition. The open loop optimal control is applied for a 
short interval of time, and feedback is introduced by reoptimization at the next time instant. This process 
presupposes that the two-point boundary-value problem can be reliably solved in a time interval that is small 
compared to the control update interval. Knowledge of the previous solution helps in providing a good starting 
point for the optimization process; however, no method has been demonstrated that can operate reliably in a 
real time environment. 

This paper examines a finite element approach to addressing this problem. Hamilton’s principle has tra- 
ditionally been used in analytical mechanics as a method of obtaining the governing equations of motion for 
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continuous systems. In the 1970’s a form of Hamilton’s principle called Hamilton’s law of varying action was 
first used by Bailey (see, for example, [2]) to obtain direct solutions to dynamics problems in the time domain, 
thus introducing Hamilton’s principle into computational mechanics. During the last decade, it was shown by 
Borri et al. [3] and by Peters and Izadpanah [4] that these direct methods, when expressed in a weak form, 
could be competitive with numerical solution of the corresponding ordinary differential equations. Later work 
by Borri et al. [5] has shown that a mixed form of the weak principle has further computational advantages, 
namely that shape functions can be chosen from a far less restrictive class of functions. 

In [6], Hodges and Bless have shown that optimal control problems can be solved in a virtually identical 
way to that of the mixed form of Hamilton’s weak principle. Hence, the method as used in [6] and the present 
paper has been called the weak Hamiltonian method for optimal control problems. Finite element methods 
have some advantages over other solution procedures; one advantage is that finite element methods provide 
the possibility for development of algorithms which converge reliably. The present method, at least for the 
problems investigated to date, is essentially self-starting. This meets a key operational requirement for on- 
board algorithms. However, application of the finite element method to optimal control problems is rather 
new. For example, Patten [7] used a Ritz-Galerkin technique with Lagragean interpolation polynomials. One 
advantage of the present formulation is the allowance for a simpler choice of shape functions. The computational 
savings which may stem from this are now under investigation. 

In this paper, a weak form governing optimal control problems is derived, and a finite element procedure is 
outlined for the solution of such problems. Numerical results for the solution of a simple trajectory optimization 
problem are presented and compared with the exact solution to demonstrate the accuracy and efficiency of 
the weak Hamiltonian finite element formulation. In anticipation of applying the present method to optimal 
guidance of a rocket booster, a simplified two-stage model suitable for this problem is presented. In [6] a 
one-stage model was analyzed and finite element results were compared to a numerical solution obtained using 
a multiple shooting method [8]. Of particular interest are the self-starting operation and the performance in 
terms of execution time and accuracy versus the number of elements used to represent the time span of the 
trajectory. 


Weak Principle for Optimal Control 

A definite analogy exists between the mixed formulation of Hamilton’s weak principle in dynamics and the 
first variation of the performance index in optimal control theory. Specifically, there is an analogy between the 
generalized coordinates and generalized momenta in dynamics and the states and co-states in optimal control 
theory. Only a brief development of the weak Hamiltonian method for optimal control problems is presented 
herein. More details on the development and the analogy with dynamics problems may be found in [6]. 


General Development 

We start with a performance index taken from Eq. (2.8.4) of Bryson and Ho [9]. Its first variation will be 
taken in a standard manner, except that states, costates, and controls will have arbitrary variations. Rather 
than setting its first variation equal to zero, however, it will be set equal to an expression which contains the 
terms that are necessary to transform all boundary conditions to the natural or “weak” type. The final weak 
form is then obtained by integration of this equation by parts in such a way that no derivatives of states or 
costates appear. 

It should be noted that the fundamental relationships are not being changed. To make certain of this, 
we will ensure that the resulting formulation produces the Euler- Lagrange equations and boundary conditions 
which have already been established in optimal control theory (see, for example, [9], Eqs. 2.8.15 - 2.8.21). 

In order to clearly understand what is meant by a “weak” formulation and the derivation of the weak 
formulation that is to follow, we first study a more simple problem. Let us start with a functional of the form 



where A and B are fixed numbers. The necessary conditions for an extremal are defined by 


( 1 ) 
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Introducing dF/di ( = / for notational convenience and integrating by parts we obtain 

J - f'^J Sydx + fSy =0 (3) 

The integrand in the above equation is the familiar Euler-Lagrange equation. The trailing term leads to the 
boundary conditions. If y is specified at x = A or x — B, then Sy — 0 at x = A or x — B respectively. This 
is referred to as a strong boundary condition. If y is not specified at one of the endpoints, then / = 0 at that 
endpoint. This is referred to as a natural or weak boundary condition. The key points to remember are that 
the trailing term itself is zero at each endpoint and that specifying y requires that 6y = 0 at a point. 

In our weak formulation, we want all the boundary conditions to be of the weak type. Thus, even if 
y(A) = 0 then Sy(A) ^ 0. To allow for this mathematically, we introduce a new variable / which represents the 
discrete value of / at an endpoint. The variation of J is now set equal to fSy |$ yielding 


6J = J dx = fSy 


This is referred to as a weak form. If we integrate by parts, we obtain 


J (jjg - f'j Sydx = (f-f)6y^ ( 5 ) 

Note that the Euler-Lagrange equation is the same as before and that the two boundary conditions are that 
/ = / at the two endpoints; thus, the trailing terms are still constrained to be zero, but in a weak sense and Sy 
need not ever vanish. 

The advantages of the weak formulation are not apparent from the above discussion. However, when we 
apply this type of formulation to problems in dynamics (see, for example [6]), or optimal control theory [6], and 
use finite elements, then we have a powerful problem-solving tool. 

Consider a system defined by a set of n states * and a set of m controls u. Furthermore, let the system 
be governed by a set of state equations of the form x = f(x, u, t). We may denote elements of the performance 
index, J, with an integrand L(x, u,t ) and a discrete function of the final states and time <j>[x(tj),tj]. In addition, 
any terminal constraints placed on the states may be placed in the set of q functions ip[x(tj), t/] and adjoined to 
the performance index by a set of q discrete Lagrange multipliers u. Finally, we will adjoin the state equations 
to the performance index with a set of Lagrange multiplier functions A(t) which are referred to as costates. This 
yields a performance index of the form 

J= f\\ T x-L-\ T J)dt-<t> -u T 1> I (6) 

Jt 0 ** 

Taking the first variation of J and setting it equal to an expression chosen so that all boundary conditions are 
of the weak type, one obtains 
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( 7 ) 


6J = [ [6\ t (x -f) + 6x T A - 6L - Sf T X] dt 

Jto 

-SxjXf - Sv T rl>\ tj 


= 6t, (A T i)| t/ - 6xl *o - 6X T (x - x) 


*/ 

to 


where 



( 8 ) 


The right hand side of Eq. (7) contains terms necessary to form all of the proper boundary conditions as natural 
ones. The quantities x and X are discrete values of the states and co-states at the initial (with subscript 0) 
and final times (with subscript /). Depending on the problem, these values will either be specified or left as 
unknowns. 

From Eq. (7), we can directly write down a weak formulation. Before this is done, however, let us examine 
this expression to ensure that it produces the correct Euler-Lagrange equations and boundary conditions. 
Integrating the 6x T term in Eq. (7) by parts and expanding the variation of L, one obtains 


-** [(£)'+ 

-“'M,-*' ( £+AT/+ ^ + *' T ^)|, ; 

"b ~ Xj^ — 6xl (*o — Aoj 

+ (*/ - */) - £Aj (* 0 - * 0 ) = 0 


( 9 ) 


where ar 0 , A 0 , */, and A j represent the values of those functions at the initial and final times, respectively. The 
coefficients of 6X T , 6x T , and 6u T in the integrand are the three correct Euler-Lagrange equations, Eqs. 2.8.15 
- 2.8.17 from [9]. There are also six trailing terms in Eq. (9) from which the boundary conditions can be 
determined. The equations corresponding to the first four and the sixth of these terms correspond to the correct 
boundary conditions in [9]. Namely, the requirement for the coefficient of 6v r to vanish yields Eq. (2.8.21). The 
requirement for the coefficient of 6t } to vanish is equivalent to Eq. (2.8.20). The requirement for the coefficient 
of 6xJ to vanish shows that the final value of A equals A/ as given in Eq.(8), which corresponds to Eq. (2.8.19). 

If A 0 is chosen as zero, the requirement for the coefficient of Sxq to vanish requires the initial value of A to equal 
zero; on the other hand, the requirement for the coefficient of 6X% to vanish requires the initial value of x to 
equal x 0 , in accordance with Eq. (2.8.18). Finally, the requirement for the coefficient of 6Xj to vanish demands 
that the final value of x equal the discrete value i/; this has no counterpart in [9] since the elements of ij are 
usually unknown. 

Having satisfied our requirement that none of the fundamental equations are altered, we may now derive 
our weak formulation from Eq. (7). In order to allow for the simplest possible shape functions when we introduce 
a finite element discretization, we do not want time derivatives of x and A to appear in the weak formulation. 
Therefore, we integrate the x term by parts in Eq. (7) yielding 
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— Sx'j X j + 8x^ Aq + 6 X'j xj — 8 A J £o — 0 


(10) 


This is the governing equation for the weak Hamiltonian method for optimal control problems. It will serve as the 
basis for the finite element discretization described below. It should be noted that normally one will encounter 
various types of inequality constraints in problems that deal with optimal control. Inequality constraints will 
be the subject of future research. 


Finite Element Solution 

Note in Eq. (10) that time derivatives of 6x and 8X are present. However, no time derivatives of x, A, u 
or 6u exist. Therefore, it is possible to implement linear shape functions for 8x and 8X within elements and 
constant shape functions for x, A, u, and within elements. 

For simplicity, let us break up the time interval into N segments of equal length At = t *Jj to • Let the values 
of time be given by U for i = 1, 2, . . . , N + 1 at the points where the time interval is broken, the so-called nodes. 
Here to = t\ and tj - tx+ Then, introduce a nondimensional elemental time r such that 


t-tj _ t-tj 
^»+i “ At 


(0<r<l) 


( 11 ) 


Now, in accordance with the above guidelines, and letting i = 1,2, . . .,N, we can choose simple linear shape 
functions 


6x = 6x(i)( 1 - r) + <fc( t+1) r 

£A = 6A(j)(l — t ) *f £A( i+1 )r 


( 12 ) 


where the arbitrary, discrete virtual states and virtual costates are defined at every node point. (The nodal 
indices are enclosed in parentheses to avoid confusion with the state column matrix index.) For the states and 
costates, we can choose even simpler shape functions 


and 


x ~ 


X = 


i 

{ 


X(i) if r = 0; 

*(,•) if 0 < r < 1; 
*(<+i) if r = 1 


A(<) if r = 0; 

A(,) if 0 < t < 1; 
A(i+i) if r = 1 


(13) 


(14) 


where in both cases, i = 1,2 For u and 6u, since their time derivatives do not appear in the formulation, 
we let 
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«=«(,) 

6u = 6&(i) 


( 15 ) 


where * = 1,2,..., TV. 

Substitution of Eqs. (12) - (15) into Eq. (10) along with Eq. (11) and 


t = U + rAt 
dt = Atdr 


(16) 


yields a set of nonlinear algebraic equations which can be assembled in accordance with standard finite element 
practice. For the sake of brevity, these rather lengthy equations for the general case are not written out explicitly 
here. They are, however, written out in [6] wherein it is noted that there are 2 n(7V + 1) equations for the states 
and costates, mN equations for the controls, q equations for the end-point constraints, and 1 equation for the 
unknown time equation. In the assembly process, £(,•) and A(,*j drop out of the equations for i = 2, ... ,7V (i.e., 
for all but the ends of the time interval to = t\ and tj = tjv+i)- Thus, the total number of unknowns is 2n(7V-|-2) 
for the states and costates, mN for the controls, q for the i/’s, and 1 for the unknown time *tv+i- Thus, there are 
2n more unknowns than there are equations, which allows for one to choose, say, x 0 in accordance with physical 
constraints and A/ in accordance with Eq. (8) and solve for the rest. The resulting equations may be explicitly 
perturbed to obtain the Jacobian and solved iteratively by a Newton- Raphson method or by any method that 
is suitable for nonlinear algebraic equations with very sparse Jacobians. 

It is also pointed out in [6] that n(7V + 1) of these algebraic equations contain the n(7V + 1) unknown 
costates and that these equations are linear in the costates. Thus, the costates can be solved for symbolically 
in terms of the other unknowns, and the remaining equations can be solved circumventing the need for initial 
estimates of the costates . This decreases the size of the problem by approximately a factor of two. 

Example: A Free-Final -Time Problem 

In this section a relatively simple optimal control problem is solved by means of the Hamiltonian weak 
formulation, and the results are compared with the exact solution. Of particular interest is the computational 
effort for varying numbers of elements, the ability of the method to converge for various values of the system 
parameters without needing new initial estimates of the unknowns, and, most important, the accuracy of the 
method versus the number of elements. 


a 



with thrust acceleration = a 


Problem Statement 

The problem is taken from [9], article 2.7, problem 9. A particle of mass m is acted upon by a thrust force 
of magnitude ma. Assuming planar motion and making use of an inertial coordinate system x x and x 2 as shown 
in Fig. 1, we may write the dynamical equations as 
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(17) 



r° 

0 

1 

°i 


f 0 


0 

0 

0 

1 

X + < 

0 

£ = 

0 

0 

0 

0 

acosu 


.0 

0 

0 

0. 


l a sin u > 


where x\ is the horizontal component of position, £2 is the vertical component of position, £3 is the horizontal 
component of velocity, and £4 is the vertical component of velocity. The control is the thrust angle ti, and 
x(0) = £0 = |_0 0 0 0J T . We want to obtain given values h of the final vertical component of position and U 

of the horizontal component of velocity. The final value of the vertical component of the velocity, must vanish, 
but we do not care what the final value of the horizontal component of position is. For the optimal control 
problem, the initial time is taken to be zero, and the final time T is to be minimized so that 0 = 0 and L = 1 
which yields J = T. The elements of the end-point constraint function 0 must vanish where 


Xf} = 


£2 — h 
x 3 -U 

£4 


1 


(18) 


In accordance with Eq. (8) 


( 0 ^ 





(19) 


Substitution of these equations into the weak form, Eq. (10), yields a system of algebraic equations whose 
size depends on N . 


Results and Discussion 

These equations are solved by a Newton-Raphson algorithm with trivial initial guesses (that are never 
changed regardless of input parameters) for N = 2. These results are then used to obtain the initial guesses for 
arbitrary N by simple interpolation. In all results obtained to date for this problem, no additional steps are 
necessary to obtain results as accurate as desired. 

Representative numerical results for £i/h versus x^/h are presented in Fig. 2 for a case with $7 = 0.75. 
Also, the control angle u versus dimensionless time t/T is presented in Fig. 3. Based on other results, not shown 
due to space limitations, accuracy for the costates is comparable to that for the states and for u. It can easily 
be seen that N = 8 gives acceptable results for all variables. 

It should be noted that the computer time on a Cyber 990 is only about 2 seconds for N = 2, N = 4, and 
N = 8 and 3 seconds for N = 16. Thus, it is relatively insensitive to N. 


Overall, the method provides very accurate results for this problem with only a few elements and for 
minimal computational effort. Furthermore, in results that are not presented herein, the Hamiltonian is seen 
to converge nicely to zero all along the trajectory as the number of elements increases. 
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Fig. 2 : Dimensionless vertical position 
x 2 /h versus horizontal position xi/h 



An Advanced Launch System Model 


In this section, a model is presented which is suitable for evaluating the potential usefulness of the weak 
Hamiltonian finite element approach in real time guidance of an advanced launch system. A two-stage vehicle 
is considered that is simplified by not allowing for any inequality constraints. 


We confine our attention to vertical plane dynamics of a vehicle flying over a spherical, non-rotating earth 
as depicted in Fig. 4. This results in the following state model for the states m (mass), h (height), E (energy 
per unit mass), and 7 (flight-path angle): 
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(20) 


m = — 


9.81/,p 

(T + qSC La \ ( V n \ 
= ( mV + 


h = V sin 7; E 


■m 


where T is the thrust, D is the drag, and V is the velocity. Here a, the angle of attack, has been adopted as a 
control variable. 



The aerodynamic, propulsion, and atmospheric models are given by the following equations: 

T =T vae - A t p 

p =p 0 (l - 0.00002255A) 5 256 for h < 11000m 
p =pn exp (- - ~ 3 ^ ) for h > 11000m 

r=R e +h ; F= y2 (i+^); 9 =^ 

( h \ v (21) 

^ =PoeXp r670oJ ; « 

D =qS [Cdo(M) + a 2 Cff a (M)] 

C La (M) =C No (M ) - C D o(M) 

a =aoVl - 0.00002255/1 for h < 11000m 
0 =295.03ms -2 for h > 11000m 

The vehicle parameters chosen for this model are based on a Saturn IB launch vehicle SA-217 [10] and are 

263.4s; I, P] = 430.4s 

Tvaci =8155800N; T vae , = 1186200N (22) 

A et = 8.47m 2 ; A t a = 5.29m 2 ; S = 33.468m 2 

where subscripts “1” and “2” refer to the first and second stages respectively. The aerodynamic coefficient data 
Coo and Cffa are presented as functions of the Mach number M in Tables 1 and 2. The physical constants 
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used in the above model are the earth’s gravitational constant p = 3.9906 x 10 14 m 3 s -2 , the earth’s mean 
radius R e = 6.378 x 10 6 m, the sea- level atmospheric pressure p 0 = 101320Nm~ 2 , the atmospheric pressure at 
11km pn = 22637Nm~ 2 , the sea-level density of air p 0 = 1.225kg m~ 3 , and the searlevel speed of sound in air 
ao = 340.3ms~ 1 . 


M 

Cd 0 

0.20 

1.00 

0.75 

0.45 

0.98* 

0.80 

1.00 

0.80 

1.02* 

0.80 

3.50 

0.20 

6.00 

0.02 


Table 1: Aerodynamic coefficient Coo versus Mach number 
(* denotes a common end point of two quadratic polynomial curves) 


M 

Cn a 

0.00 

6.20 

0.50 

6.35 

0.98* 

7.70 

1.00 

7.70 

1.02* 

7.70 

2.50 

5.20 

4.40* 

4.70 

5.00 

5.50 

6.00 

6.00 


Table 2: Aerodynamic coefficient C^a versus Mach number 
(* denotes a common end point of two quadratic polynomial curves) 


The performance index is 


J = = m| (/ 


(23) 


and the final time tj is open. The initial conditions specified are m(0) = 5.2 x 10 5 kg, h(0) = 1800m, E( 0) = 
—6.25 x 10 7 m 2 s~ 2 , and 7(0) = 75°. The final energy is specified as E{tj) = —1.25 x 10 7 m 2 s -2 . The burnout 
mass of the first stage is 192000 kg and the drop-mass of the booster is 51000 kg. 

For this two-stage model, we must modify our formulation somewhat. We must accomodate for the unknown 
staging time t,, the constraint on the mass at t, (as opposed to a constraint on the states at the final time), 
the jump in the mass at t,, the jump in the mass costate at t,, the condition for a continuous Hamiltonian at 
t, (continuous since <f> and ip are not explicit functions of time), and finally the change in state equations at 
t, (due to the change in thrust). We further point out that the control u is discontinuous at the staging time. 
However, since only discrete mid-point values of u are solved for, the jumps are allowed to occur automatically 
at the nodes. 

The new performance index (with L — 0) is 


= [A T (* - fi) & + J* [A T (* - h) dt - vi V>i 


-<t>\ 




t. 


(24) 


351 



where f\ and / 2 represent the state equations before and after t s respectively, ) — 192000 and 

= E(tf)+ 12500000. 

The weak formulation is derived exactly as before and the same shape functions can be employed. This 
time, however, we will discretize the time from 0 to t s with N\ elements and the time from t s to tj with N 2 
elements. The algebraic equations shown in [6] for the weak formulation are readily modified to account for the 
various jumps and state equation discontinuities. The resulting equations are 


A i +SXJ 


“ — (/i)< 


-SuT |a<i (H) T ^] }-«. [A(ti) - H(t:)} - 6* (rh- Ni+i - 192000) 

+6x\\q — SXfxo — + l + l ~ ^JVj + l) + ^JVi + l (*AT, + 1 ~ *JVi + i) 


At 2 li 

2 


(/>), 


[a. - & (f )' a] - «&, [* + t* CM. 

-Mf [a.j (f) aJ J-<1; [«(</)] - in (e, + 12500000) 
+^A^ 1+A r 3+1 (i/) — ) = 0 


(25) 


From Eqs. (8) and (23) it is seen that \j is given by [1 0 v 2 0j T . Also, we note that the only jumps 

are in the mass state and the mass costate and these jumps are 


m(t;)-m(t+) = 51000 

(C 5 ) — A m (<+) = v \ 


The finite element equations are solved using the method of Levenberg-Marquardt as coded in the IMSL 
subroutine ZXSSQ [11]. Running a case for a few elements generates a good approximation for larger numbers 
of elements. Initial guesses do not need to be very accurate, but the method is not nearly as computationally 
efficient as a Newton- Raphson procedure where sparsity in the Jacobian could be exploited. 


In Figs. 5 and 6 numerical results for the ALS model are given for 4 and 8 elements in each time interval. 
(The number of elements in each interval is completely arbitrary.) In Fig. 5 the altitude profile is shown and 
the control history is shown in Fig. 6. From past experience with the one-stage model, we believe the N\ = 8 
and N 2 = 8 to be a converged result. Furthermore, we see that even the Ni = 4 and N 2 = A result gives a 
reasonable approximation to the solution. It should be noted that these results are not realistic because of the 
absence of state constraints, and because we have large angles of attack (more than 30° at some points) even 
though we assumed small singles in the state equations. However, they do suffice to illustrate the power of the 
method. 


352 




As an indication of the accuracy of the method in a global sense, the Hamiltonian was observed to converge 
to zero (the exact answer) all along the trajectory. The finite element results are converging to the exact solution 
as N increases. 


Conclusion 

In this paper we present a weak Hamiltonian formulation for optimal control problems. Results are pre- 
sented based on the weak Hamiltonian finite element formulation for a simple optimal control problem which 
show the method is reliable, efficient, and accurate. For this and other problems it is easily programmed with 
a self-starting algorithm. 

To address the future needs of real-time guidance for future space transportation systems, we present a 
four-state model for the dynamics of mass, altitude, energy, and flight-path angle. The angle of attack is 
the control. The results show the power and efficacy of the present approach. It has several advantages for 
applications in real-time guidance. It is not only reliable and efficient but has excellent self-starting capability. 
Furthermore, initial guesses of the costates are not needed, and the method exhibits a graceful degradation in 
performance with reduction in number of elements. 

For future research we intend to concentrate on improving computational efficiency by exploitation of the 
relatively significant level of sparsity in the Jacobian. We also will incorporate inequality constraints, and begin 
to address the avoidance of atmospheric anomalies. 
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